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Abstract 


Nonlinear  Scale-based  signal  and  analysis  techniques  are  widely  viewed  as  a  viable  alter¬ 
native  to  meet  the  challenges  of  increasingly  complex  engineering  problems.  The  research 
effort  described  in  this  report  addresses  two  problems  which  require  such  techniques: 

1.  To  make  Global  Positioning  System  more  robust  to  outside  interference  and  pertur¬ 
bations 

2.  To  perform  filtering  on  images  for  Automatic  Target  Recognition  driven  by  features 
which  may  be  potential  signatures. 


Chapter  1 

Summary  of  contributions 


The  funding  of  AFOSR-F49620-98-1-0190  grant  has  helped  the  PI  establish  a  solid  footing 
for  a  new  research  program,  and  has  resulted  in  su bst anti aftcont r i bu t ions  in  Nonlinear 
Multi-scale  Estimation  and  Filtering.  A  Master’s  thesis  funded  at  100%  authored  by  B. 
Karacali  was  completed  [1]  and  two  partially  (50%  funding  each)  supported  Ph.D.  theses 
by  Y.  F.  Bao  (May’01)  and  O.V.  Poliannikov  (May  02)  are  near  completion  .  Partial 
support  of  this  grant  has  also  led  to  the  publication  of  four  journal  papers  [2,  3,  4,  5]  and 
three  journal  papers  under  review  [6,  7,  8],  as  well  as  two  invited  book  chapters[9,  10]  and 
several  conference  papers  at  ICASSP  and  ICIP  [11,  12,  13,  14]  and  invited  national  and 
international  seminars,  including  an  invitation  for  a  one  month  visit  in  dhe  department  of 
information  sciences  at  University  of  Tokyo  (Japan)  (Nov.  1999),  courtesy  of  the  Japanese 
Society  for  Advancement  of  Science,  and  an  invitation  as  principal  guest  speaker  at  a 
workshop  on  Mathematics  in  Imaging  at  University  of  Hong  Kong  (China)  (Dec.  2000). 

In  the  course  of  this  project,  working  contacts  have  also  been  initiated  with  key  technical 
members  of  the  GPS  group  of  the  Air  Force  Laboratories.  Our  plan  in  the  very  near 
future  is  to  contact  the  ATR  groups  around  the  Air  Force  Laboratories  and  to  better 
communicate  some  of  our  results  and  help  fine  tune  them  to  specific  problems  of  direct 
interest  to  the  Air  Force. 

The  increasing  complexity  of  engineering  problems  in  general,  together  with  the  demand 
for  additional  performance  in  signal  and  image  processing  problems  in  particular,  have  led 
to  an  intensification  of  search  for  nonlinear  alternative  approaches.  This  is  widely  viewed 
as  a  promising  framework  which  would  provide  reliable  as  well  as  robust  techniques  to 
meet  new  challenges. 

Towards  that  goal,  the  work  performed  over  the  duration  of  the  current  grant  has  focused 
on  multi-scale  nonlinear  techniques  with  two  main  application  thrusts: 

1.  Robust  methods  for  Global  Positioning  System  applications 

2.  Automatic  Target  Recognition  in  Synthetic  Aperture  Radar,  Infra-red  and  Inverse 
Synthetic  Aperture  Radar  imaging. 


In  GPS.  a  frequently  encountered  important  problem  in  particularly  adverse  conditions 
(e.g.  intentional  jamming)  is  a  loss  of  lock  which  in  turn,  may  lead  to  unexpected  delays 
and  hence  to  serious  consequences  in  critical  GPS  applications. 

The  first  part  of  this  report  addresses  such  a  problem,  using  nonlinear  wavelet  techniques 
which  are  well  known  to  show  particular  resilience  to  noisy  environments  and  and  to 
exhibit  robustness  to  perturbations.  Specifically,  we  develop  a  new  code  lock  detection 
technique  which  among  others,  facilitates  the  detection  of  the  location  of  the  data  bit 
transition.  In  the  second  part  of  this  report,  we  explore  the  potential  of  nonlinear  scale- 
based  techniques  in  ATR.  Specifically,  we  develop  a  probabilistic  framework  for  well  known 
nonlinear  filtering/diffusion  techniques  which  give  an  intuitively  appealing  interpretation 
which,  we  in  turn  use  to  develop  more  performing  techniques  in  blind  image  enhancement 
[6]  and  techniques  which  use  features  inherent  to  the  image  to  drive  the  filtering. 


Part  I 
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Robustness  of  GPS 
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Chapter  2 


Nonlinear  Wavelet  Techniques  in 
GPS 


2.1  INTRODUCTION 

Global  Positioning  System  (GPS)  is  today’s  most  advanced  navigation  system,  which 
incorporates  a  position  determination  tool  through  a  web  of  satellites  that  broadcast 
their  exact,  position  in  orbit  repetitively  at  particular  time  instants.  By  decoding  each 
satellite’s  message  to  get  their  exact  position,  and  by  computing  their  individual  distances 
to  the  receiver,  a  user  can  determine  their  position  on  the  earth’s  surface  with  very  good 
accuracy. 

A  successful  operation  of  the  receiver  consists  of  its  simultaneous  decoding  of  at  least  four 
different  satellites  which  in  turn  requires  tracking  of  the  signal  carrier  frequency  and  code 
delay.  This  requirement,  as  illustrated  in  Figure  2.1,  is  reflected  by  the  demodulation 
step  as  well  as  the  receiver’s  alignment  of  its  Pseudo  Random  Noise  (PRN)  sequence  to 
that  of  the  received  signal.  On  account  of  currently  used  discriminators’  vulnerability  to 
bit  transitions,  the  receiver  has  to  acquire  lock  in  both  domains  within  a  20  m-sec  time 
interval. 

The  present  work  addresses  this  problem  and  provides  additional  robustness  to  carrier 
frequency  and  code  delay  lock  discriminators  which  are  subject  to  such  navigation  message 
bit  transitions.  In  what  follows,  we  first  summarize  the  GPS  receiver  operation,  and 
subsequently  describe  the  robust  discriminators  together  with  performance  comparisons 
with  the  current  discriminators  using  real  GPS  data. 
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2.2  GPS  OVERVIEW 


Figure  2.1  shows  the  block  diagram  of  a  generic  GPS  receiver.  Upon  its  preconditioning 
and  discretization  at  a  lower  intermediate  frequency  (IF),  the  radio  frequency  (RF)  signal 
is  fed  to  receiver  channels  which  track  signals  from  different  satellites. 


#1 


Figure  2.1:  The  block  diagram  of  a  generic  GPS  receiver 

The  receiver  channels  provide  the  demodulated  navigation  messages  of  all  observable 
satellites  to  the  receiver  processing,  which  then  feeds  the  readings  to  navigation  processing 
that  computes  the  position  for  the  user. 

The  accurate  computation  of  the  position  is  therefore  crucially  dependent  on  accurate 
demodulation  of  several  satellites’  signals,  which  is  intrinsic  to  the  receiver  channels. 


2.2.1  Generic  Receiver  Channel 

A  simplified  block  diagram  of  a  generic  receiver  channel  is  shown  in  Figure  2.2.  The 
incoming  digital  IF  signal  is  first  stripped  off  from  its  carrier  through  the  sine  and  cosine 
maps,  and  then  collapsed  to  baseband  once  correlated  with  the  receiver  generated  PRN 
code. 

The  receiver  processor  acquires  frequency,  phase  and  code  delay  errors  through  Ie,  Ip,  I;, 
Qe,  Qp  and  Q/  readings,  and  tries  to  preserve  the  locks  by  adjusting  the  frequency  and 
phase  of  the  sine  and  cosine  maps,  and  the  clock  rate  of  the  PRN  code  generator. 
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2.3  CURRENT  LOCK  DISCRIMINATORS 

The  lock  discriminators  used  in  this  work  are  summarized  below.  A  more  through  pre¬ 
sentation  is  available  in  [15]  §5.1. 

2.3.1  Frequency  and  Phase  Lock  Discriminators 

In  order  to  be  able  to  strip  the  carrier  from  the  spread  spectrum  signal,  the  receiver  needs 
to  know  the  precise  carrier  frequency  of  the  incoming  digital  IF  signal.  One  of  the  most 
commonly  used  frequency  lock  discriminators  is  given  by: 

_  -  02 

^  (t2  —  *i)2n 

where  0*  =  atan2(I*, Q*),  and  I*  is  the  lp  sample  observed  at  time  U,  for  i  =  1,2.  This 
discriminator  gives  the  true  frequency  error  in  a  neighborhood  inversely  proportional  to 
the  integration  time.  The  plot  of  the  discriminator  output  versus  true  error  is  shown  in 
Figure  2.3.  The  dashed  line  in  Figure  2.3  indicates  the  ideal  discriminator  output. 


Figure  2.3:  The  error  output  of  the  regular  FLL  discriminator 

The  phase  lock  is  generally  obtained  by  a  Costas  loop.  The  optimal  phase  lock  discrimi¬ 
nator  [15]  is  given  by 

5<p  =  atan(Ip/Qp). 
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This  phase  discriminator  gives  the  true  phase  error  in  a  neighborhood  of  ±11/2. 


j* 


2.3.2  Code  Delay  Lock  Discriminators 

Among  a  number  of  common  delay  lock  discriminators,  the  normalized  early  minus  late 
envelope  discriminator  given  by 


E  \/ifTQ!  +  £  \/SFTq? 

provides  an  amplitude  insensitive  delay  error  within  ±  0.5  bit  delay.  In  spite  of  its  high 
computational  load,  it  will  be  used  in  this  text  as  a  basis  for  comparison  for  the  proposed 
code  delay  lock  discriminator.  The  theoretical  error  output  of  this  common  discriminator 
is  showm  in  Figure  2.4. 


Figure  2.4:  The  theoretical  error  output  of  the  normalized  power  DLL  discriminator 


2.4  PROPOSED  LOCK  DISCRIMINATORS 

The  lock  discriminators  described  here  are  based  on  research  work  presented  in  [16].  In 
[16]  §5.2.2,  an  operator  Oj0  acting  on  a  signal  s  that  returns  high  values  for  signal  regions 
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with  embedded  singularities  is  defined.  An  operator  d  defined  through  Oj0  with  jo  =  1  as 

d(s)  =  max  |max{Oi(s)},  max{Oi(s)}| 

where  k  is  the  time  index  of  Qj0(-)\  s  :  [0,  to]  — >•  K,  and  s(t)  is  given  by 

gU\  _  /  s(t)  ift  <=  [°>  |] 

'  '  [  —  s(t)  otherwise  ’ 

is  shown  to  be  robust  to  sign  transitions  in  s(t),  which  creates  potential  for  use  in  lock 
discriminators  with  robustness  to  bit  transitions. 


2.4.1  Robust  Frequency  Lock  Discriminator 

Let  Ip  denote  the  input  of  the  integrate  and  dump  operator  that  produces  Ip  samples. 
Then,  the  robust  frequency  loop  lock  (FLL)  discriminator  is_defined  as 

Tf  <Pl-<P2 

Jr  (h-h)  2n 

where  ipi  =  atan2(/*,  Qlp)  with  I‘lp  =  d(lp)  ■  sign(]P  lp)  observed  for  a  certain  period  of 
time  until  U,  for  %  =  1,  2.  Figure  2.5  shows  the  discriminator  output  versus  true  frequency 
error  to  be  almost  linear  with  slope  1  around  the  true  frequency,  indicated  by  the  dashed 
line,  and  therefore,  the  discriminator  indeed  can  be  used  to  determine  the  frequency  error 
between  the  incoming  and  the  receiver  generated  carriers. 


2.4.2  Robust  Code  Delay  Lock  Discriminator 

The  robust  code  delay  lock  (CDL)  discriminator  is  defined  by  the  following  equation: 

err  -  V^(Se)  ~  V^(s') 

\J  d(se)  +  a/  d(si) 

where 

s  f  ydf~+~Qf  •  sign(Ie)  if Ej>Eq 
\  Qe  '  sign(Qe)  otherwise 

Eq  =  EQI- 

The  empirical  plot  of  the  error  output  versus  input  delay  of  the  robust  discriminator  is 
shown  in  Figure  2.6.  Note  that  the  plot  is  very  similar  to  the  theoretical  input  output 
response  of  the  regular  code  delay  lock  detector,  showing  linear  error  within  ±  0.5  bit  of 
delay  error. 

1Note  that  k  =  1, . . .  ,  2-'°. 
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2.5  PERFORMANCE  ANALYSIS 


The  performance  analysis  of  the  regular  lock  discriminators  and  the  proposed  robust 
discriminators  is  conducted  using  real  GPS  data2.  The  data  consists  of  the  output  of  the 
analog-to-digital  converter  in  a  GPS  receiver,  which  is  the  input  fed  to  individual  receiver 
channels  for  carrier  and  code  stripping  processes. 

2.5.1  Frequency  Tracking  Performance 

The  vulnerability  of  the  regular  frequency  lock  loop  discriminators  to  bit  transitions  in 
the  navigation  signal  is  well  known.  This  is  also  demonstrated  in  Figure  2.7  which  shows 
frequencies  at  which  the  regular  discriminator  indicates  the  frequency  lock.  For  some  time 
where  no  bit  transitions  occur  in  the  data  stream,  the  lock  appears  stable.  Any  inadvertent 
bit  transition,  however,  causes  the  discriminator  to  either  lock  at  another  frequency,  or 
loose  the  frequency  lock  altogether.  Note  that  the  spikes  \^nch  depict  the  loss  of  the 
frequency  lock  are  exactly  20  msec  apart,  indicating  that  they  indeed  correspond  to  the 
bit  transitions  of  the  50  Hz  navigation  signal. 

The  proposed  robust  FLL  discriminator  however,  does  not  show  erratic  behaviors  around 
the  bit  transition  regions  that  fool  the  regular  discriminator,  maintaining  the  frequency 
lock  regardless  of  the  navigation  message  bit  transitions,  as  shown  with  solid  lines  in  Figure 
2.8.  The  dashed  lines  in  Figure  2.8  represent  the  regular  frequency  lock  discriminator 
track. 


2.5.2  Code  Delay  Tracking  Performance 

In  [16],  it  was  empirically  shown  that  the  robustness  to  bit  transitions  in  the  conventional 
code  delay  lock  detectors  based  on  envelope  detection  was  at  a  cost  of  a  higher  sensitivity 
to  noise.  The  performance  comparison  between  the  regular  code  delay  lock  discriminators 
and  the  proposed  discriminators  should  hence  be  carried  out  as  a  function  of  prevailing 
noise. 

Towards  that  end,  the  real  GPS  data  is  demodulated  separately  using  regular  and  pro¬ 
posed  delay  lock  discriminators,  both  with  the  regular  frequency  lock  discriminator,  on  a 
signal  region  with  no  navigation  message  bit  transitions.  This  procedure  is  fifty  hundred 
times,  and  with  various  levels  of  white  noise  added  artificially  to  the  initial  real  GPS 
signal.  The  results  are  shown  in  Table  2.1. 

The  simulation  results  indicate  that  in  similar  conditions,  the  proposed  delay  lock  dis¬ 
criminator  offers  better  noise  performance  than  the  regular  discriminator,  under  additive 
white  Gaussian  noise.  The  increasing  noise  levels  impair  both  discriminators,  but  at  SNR 

2We  would  like  to  thank  Michael  S.  Braasch  and  Maarten  U.  de  Haag  of  University  of  Ohio  for 
providing  the  real  GPS  data. 
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frequent 


Figure  2.7:  The  regular  FLL  discriminator’s  frequency  tracking 


Table  2.1: 


SNR 


20  15  10  5 


0 


> 


Figure  2.8:  The  robust  FLL  discriminator’s  frequency  tracking 


levels  where  the  regular  discriminator  looses  the  signal,  the  robust  discriminator  still  offers 
accurate  signal  tracking. 


2.6  CONCLUSION 

So  far,  we  have  presented  a  performance  analysis  on  the  proposed  frequency  lock  and  code 
delay  lock  discriminators  on  real  life  GPS  data,  and  compared  their  performances  to  two 
existing  discriminators.  The  proposed  frequency  lock  discriminator  proved  to  be  robust 
to  navigation  data  bit  transitions.  This  is  a  big  problem  for  regular  discriminators  as  it 
limits  the  predetection  integration  time  and  the  noise  performance  thereafter.  The  use  of 
a  robust  frequency  lock  discriminator  therefore  should  alleviate  the  20  msec  restriction  on 
the  predetection  integration  time,  making  operation  possible  in  much  heavier  attenuation, 
noise  or  interference  environments. 

In  comparison  to  envelope  detection,  the  proposed  code  delay  lock  discriminator  showed 
better  performance  in  establishing  a  delay  lock  under  real  noise  conditions  This  confirms 
the  simulated  results  presented  in  [16],  and  offers  promising  functionality  of  the  newly 
proposed  techniques  in  adverse  signal  conditions. 
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Part  II 

* 

Nonlinear  Diffusion:A  Probabilistic 

Approach 
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Chapter  3 

Anisotropic  Diffusion 


3.1  Introduction 


T 


Scale-based  analysis  has  played  an  increasingly  important  role  in  signal  and  image  anal¬ 
ysis  since  Witkin’s  ground  breaking  paper[17],  in  which  a  so-called  linear  scale  space  was 
constructed.  This  was  based  on  the  conclusion  that  convolving  a  signal  with  a  Gaussian 
kernel  was  equivalent  to  evolving  it  with  a  Heat  differential  operator  where  time  plays  the 
role  of  scale  [18,  19].  This  linear  filtering  approach,  however,  presented  a  major  limita¬ 
tion,  namely  that  important  details  in  signal/image  also  get  smoothed  away  along  with 
the  noise.  On  the  other  hand,  and  with  a  different  twist,  Mallat  [20]  proposed  a  system¬ 
atic  nonlinear  multi-scale  analysis  framework  using  wavelet  bases.  The  computational 
efficiency  of  a  wavelet  analysis  together  with  its  ability  to  focus  on  and  localize  salient 
features  of  a  signal  via  its  multi-scale  coefficients  have  secured  it  great  popularity.  This 
framework  provided  a  solid  foundation  for  a  number  of  subsequent  studies  on  De-noising, 
segmentation,  etc.  [21,  22].  The  implicit/explicit  assumption  of  independence  of  the 
wavelet  coefficients  in  the  analyses  of  these  De-noising  techniques  and  others  turned  out' 
to  be  a  limitation  in  some  cases  and  degraded  their  performance  and  overall  effectiveness 
in  others.  An  attempt  to  mitigate  such  a  problem  was  provided  for  1-dimensional  signals 
in  [23]  by  accounting  for  some  Markovian  structure,  the  solution  appeared  to  improve 
performance  albeit  at  a  substantial  computational  cost.  Its  applicability  for  images,  to 
the  authors’  best  knowledge,  was  never  given.  Given  the  intrinsic  ability  of  the  continuous 
scale  approach  to  preserve  intra-scale  correlation  information  as  well  as  inter-scale  infor¬ 
mation,  it  is  reasonable  to  expect  techniques  based  on  such  approach  to  result  in  more 
efficient  and  effective  filtering.  Such  a  clever  notion  was  first  used  by  Perona  and  Malik 
in  their  landmark  paper  [24]  where  they  aimed  at  preserving  important  sharp  features 
such  as  edges  in  images.  Their  technique  may  also  be  viewed  as  a  nonlinear  filter  whose 
selective  smoothing  is  based  upon  the  computed  local  gradient  (maximal  smoothing  in  low 
gradient  or  homogeneous  regions,  and  minimal  smoothing  in  high  gradient  regions).  The 
novelty  of  this  approach  together  with  its  most  promising  results  triggered  a  tremendous 
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research  activity  in  computer  vision  and  applied  mathematics  (see  [25]  for  a  good  review 
of  the  literature  as  well  some  other  fundamental  papers  such  as  [26,  27]). 

A  number  of  very  good  papers  have  recently  provided  inspiring  variational  interpretations 
to  various  nonlinear  smoothing  techniques  [25].  Specifically,  many  existing  nonlinear 
evolution  equations  were  shown  to  result  from  a  minimization  of  energy  functionals  whose 
Euler-Lagrange  equations  led  to  a  gradient  descent  method,  thereby  giving  rise  to  a  PDE. 
All  of  these  evolution  techniques  required  an  a  priori  stopping  time,  the  absence  of  which 
almost  always  led  to,  as  is  well  known,  to  a  complete  smoothing  of  the  signal/image 
(i.e.,  the  steady  state  of  the  PDE).  Furthermore  most  of  the  existing  developments,  if 
not  all,  have  been  predominantly  deterministic  in  nature,  with  little  to  no  stochastic 
treatment  or  interpretation  of  the  diffusion.  The  characteristics  of  the  random  process 
which  underlies  a  diffusion  have  hence  been  overlooked,  and  their  overall  influence  on  the 
solution  in  different  scenarios  has  remained  unclear.  Poliak  et.  al.  [5]  recently  proposed 
an  approach  addressing  robustness  issues,  and  showed  some  remarkable  results  for  a  wide 
class  of  perturbation  noises.  The  analysis  remained  as  in  all  other  cases,  fundamentally 
deterministic,  and  also  required  knowledge  of  the  stopping  ti&e  for  the  evolution. 

Towards  understanding  the  intrinsic  stochastic  behavior  of  nonlinear  diffusions,  we  adopt 
a  probabilistic  view  of  diffusion  (and  partially  described)  in  [11],  and  propose  a  framework 
where  nonlinear  diffusions  are  cast  and  are  given  an  insightful  interpretation.  Additional 
extensions  are  also  proposed  in  [12].  In  this  report  we  describe  in  sufficient  detail  an 
alternative  view  of  nonlinear  filtering/diffusion  which  may  also  be  found  in  [6].  This  in 
fact  is  instrumental  in  our  providing  an  alternative  interpretation  opexisting  methods, 
e.g.,  Perona-Malik  equation,  and  in  using  the  gained  insight  to  propose  a  solution  to  its 
well  known  limitations.  More  specifically,  we  view  an  evolution  equation  as  a  solution  to 
a  controlled  diffusion  [28]  resulting  from  an  optimization  of  an  energy  functional.  This 
ultimately  leads  to  a  two  to  four  state  Markov  Chain  (MC)  with  one  step  transition 
probabilities  well  adapted  to  preserving  the  salient  features  of  an  image/signal  (  such 
as  edges)  while  smoothing  away  the  noise.  As  will  be  elaborated  on  further  below,  in 
addition  to  a  marked  performance  improvement  over  P-M  equation,  and  by  way  of  our 
newly  proposed  technique  we  are  able  to  lift  a  longstanding  problem  in  nonlinear  diffusion , 
namely  requiring  prior  knowledge  of  a  stopping  time.  We  in  fact  show  that  the  stable  point 
for  our  equation  is  a  staircase  function.  The  immediate  applications  of  such  a  filtering 
technique  are  signal/image  enhancement  as  well  as  segmentation  and  feature  extraction 

[7]- 

Some  background  material  is  first  presented  in  the  next  section  and  prior  to  our  reformu¬ 
lation  of  the  diffusion  problem  in  Section  3.  Following  a  discrete  formulation,  we  proceed 
to  alternatively  reinterpret  in  Section  4  nonlinear  diffusion  equations  as  non-homogeneous 
controlled  Markov  Chains.  In  Section  5,  we  use  the  insight  gained  in  the  previous  sections 
to  propose  a  new  algorithm  whose  evaluation  through  numerous  substantiating  examples 
is  provided  in  Section  6.  We  finally  give  some  concluding  remarks  in  Section  7. 
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3.2  Probabilistic  View  of  Diffusion 


A  stochastic  process  Xt  may  be  defined  as  a  parameterized  collection  of  random  variables 
{Xt}teT  defined  on  a  probability  space  (Q.  T ,  P)  and  assuming  values  in  Kn  in  general, 
so  that, 


Vcu  e  Q,  u>  — ^  A  t  £  T , 

where  Q  is  the  usual  sample  space,  T  the  a-field  and  P  the  probability  measure.  A  nice 
and  intuitively  appealing  interpretation  for  u>  is  that  of  a  particle  whose  position  at  time 
t  is  given  by  Xt(u). 

Definition  1.  A  stochastic  process  which  satisfies  the  following  stochastic  differential 
equation(SDE)  is  called  an  Ito-diffusion[29j. 

dXt  =  b(t,  Xt)dt  +  aft,  Xt)dBu  (3.1) 

where  Bt  is  a  standard  Brownian  motion  and  b(t,x),a(t,  x)  (fire  drift  and  diffusion  coeffi¬ 
cients,  (d  =  1, 2  in  this  paper).  Furthermore,  if  the  drift  and  diffusion  terms  b(t,  x),a(t ,  x) 
are  associated  with  a  function  v(t,x),  Xf  =  Xt  is  called  a  controlled  diffusion,  where  Xf 
satisfies  the  SDE 


dXvt  =  b{t,  v{t,  Xt))dt  +  a(t,  v{t,  Xt))dBt. 


(3.2) 


Denote  by  p(s,x,t,dy )  the  probability  transition  function  of  a  stochastic  process  Xt,  i.e., 
p(s,x,t,dy)  =  P{Xt  e  dy\Xs  =  x),  and  by  p(s,x,t,y )  the  probability  transition  density 
function  of  a  particle  starting  at  location  x  and  time  s  and  reaching  y  at  time  t.  The 
transition  probability  is  normally  determined  by  the  drift  and  diffusion  coefficients,  which 
characterize  how  the  diffusion  behaves.  If  there  is  a  control  factor  in  the  coefficients, 
it  is  reflected  by  a  constrained  transition  probability  for  the  particle.  The  infinitesimal 
generator  (i.e.,  continuous  operator  which  described  such  a  motion)  of  the  diffusion  in 
Eq.  (3.2)  can  be  then  written  as  : 

ij=l  1  3  1 


where  a (t,v(t,x))  =  a(t,v(t,x))Ta(t,v(t,x)),  The  solution  of  the  resulting  PDE 


dUt{x) 

dt 

U(s,  x) 


LvtUt(x) 

f(x)  for  some  s  >  0, 


(3.3) 


can  be  written  as  an  expected  value  Ut(x)  =  E x{f(Xt)}  =  f  f(y)P(s,x,t,dy)  [30].  If  we 
let  s  =  0.,  Eq.(  3.3)  is  the  so-called  Markov  forward  equation.  A  backward  equation  can 
be  viewed  as  a  reverse  time  process,  with  the  following  form: 
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Figure  3.1:  A  particle  (pixel)  may  diffuse  over  many  possible  paths,  and  an  average  is 
usually  computed. 


+LTU,(x)  =0 

Ut(x)  —  f{x)  for  some  T  >  0,  (3.4) 

where  Lvt*(-)  is  an  adjoint  operator  of  Lvt.  The  solution  can  again  be  expressed  as  Ut(x)  = 
Et>x(f  (Xt)) ,  where  T  is  the  fixed  terminal/end  time  of  the  diffusion  (or  initial  in  the 
case  of  an  inverse  diffusion).  Note  that  a  backward  diffusion  equation  can  be  viewed  as  a 
forward  diffusion  by  merely  selecting  t'  =  T  —  t. 

An  illustrating  example  of  a  non-controlled  diffusion  is  the  process  described  by  the  PDE 
in  Eq.(  3.8),  in  which  Lvt  is  specified  as  a  Laplacian  operator  A  (i.e.,  ^  +  J^-).  Diffusion 
of  heat  in  a  homogeneous  medium  fundamentally  stems  from  the  motion  of  particles.  It 
can  be  shown  that  the  inherent  randomness  of  this  motion  is  well-described  by  a  Brownian 
motion  Bt  [29],  where  an  individual  outcome  to  £  D  in  the  prevailing  sample  space,  may  be 
associated  to  a  particle.  The  process  Bt  can  then  be  interpreted  as  the  distance  traveled  by 
particle  u>  at  time  t.  It  is  well  known  that  such  a  transition  density  for^a  Brownian  motion 

in  1-D  case,  for  instance,  is  a  Gaussian  PDF  p(t,x,y)  =  r  Vrr,  y  e7l,t>0 

(recall  that  a  Brownian  motion  has  independent  Gaussian  increments).  It  is  thus  clear 
that  a  stochastic  interpretation  of  a  solution  (if  it  exists)  subject  to  some  differentiability 
conditions,  can  be  given  by  way  of  an  ensemble  average  [30] 

Ut(x)  =  Ex{f(Bt)}  =  [  p{t,x,y)f(y)dy,  (3.5) 

JR 

where  the  expectation  E(-)  is  computed  over  all  possible  reachable  positions  x\  starting 
at  position  x.  In  the  2-D  case,  it  is  similarly  possible  to  have  such  an  interpretation  as 
displayed  in  Fig.  3.1.  The  times  t  —  fj  and  t  =  t2  are  the  instants  at  which  all  possible 
positions  are  averaged  to  yield  a  solution  at  the  respective  times. 


3.3  Diffusion  on  a  Lattice 

As  previously  noted,  our  chief  interest  here  is  to  propose  a  framework  within  which  a 
stochastic  interpretation  of  a  diffusion  (or  more  generally  of  the  so-called  scale  space 
analysis)  is  achieved,  and  is  in  turn  instrumental  in  gaining  insight.  Towards  that  end 
and  to  further  extend  and  possibly  improve  on  existing  techniques,  we  begin  by  discretizing 
the  space  as  well  as  the  scale/time  variables. 
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3.3.1  Discrete  Approximation  of  Diffusion 

Recall  that  a  symmetric  one-dimensional  (1-D)  random  walk  is  well  known  to  converge 
to  a  Brownian  motion  as  r  — *  0  and  5—^0,  with  t,  8  respectively  denoting  scale  and 
distance  discrete  steps.  A  particle  following  such  a  trajectory  will  move  on  a  1-D  lattice 
with  probability  1/2  to  the  left  or  to  the  right,  while  on  a  2-D  plane,  it  will  move  to 
any  of  the  four  nearest  neighbors  (east,  west,  north,  south)  with  equal  probability  of  1/4. 
Formally,  in  2-D  space,  we  write  the  spatial  variable  (a; j ,  x\)  =  (ay  -H<5,  x2+i8)  with  i  e  Z 
and  the  scale  tn  =  nr  with  n  G  N,  and  we  denote  the  one  step  transition  probability  of  a 
particle  from  (a;?, a:®)  to  (a^aq)  at  the  nth  scale  step,  by  pn((a;°,a;2),  (xi,x2)).  As  a  result, 
we  obtain  a  standard  form  from  Eq.  (3.5),  namely  the  probability  of  a  particle  being  at 
(aq,  x2)  at  scale/time  (n  +  l)st  step  as  r  — »  0  and  8  — >  0, 

Proposition  1.  The  following  discrete  equation, 

pn+i((a^i,^2))  (ah,  £2))  =  \Pn({xi,x°7),{x-8,y)  +  \pn((x0l,x02),{xi+8,x2)) 

+4Pn((*i>*2)>  (*i»®2  -  <5))  +  (xi,x2  +  5)) 

(3-6) 


converges  to 


dptifA^ 2),  OP,  3-2)) 

dt 


a2pt((xg,a:^),  (aq,^)) 

dx\ 


dx\ 


(3.7) 


Proof:  Subtracting  p^x^x®)  from  both  sides  of  Eq.  (3.6)  and  dividing  it  by  r,  we  obtain 


[pn+i((x$,a;^),  (aq,aq.))  -Pn({x°i,xl),  (xux2 ))]/r  = 
if  bn((®i,*5),  (®i  -  (S,^))  -  2pn((x°1,x°2),  (aq,aq))  +  p„((zS,a:!j),  (a*  +  <5,ar2))]  + 
if  [Pn((®l,®a)i  (®1,®2  -  <5))  -  ^Pn({Xi,X°2),  (®  1,0:2))  + 

Pn({Xi,X2)>  {X1,X2  +  5))] 


which  upon  letting  r  =  <52/ 4  and  <5  ->■  0,  concludes  the  proof.  ■ 

The  numerical  approximation  of  a  linear  diffusion  is  hence  readily  implemented  by  way 
of  a  random  walk. 

When  considering  a  1-D  Brownian  motion  on  a  compact  interval,  a  reflecting  wall  should 
be  accounted  for,  making  the  Markov  Chain  (MC)  aperiodic  and  recurrent  and  for  which 
a  solution  to  A U(x)  =  0 ,Uq(x)  =  f(x)  takes  the  form  U(x)  =  Eg(f( XT)).  This  also 
implies  a  discrete  solution  U(x)  =  '^Tpif(Xi)  where  Pi  is  the  probability  of  a  particle  to 

i 

be  in  state  i  as  t  grows  large.  The  independence  of  the  solution  U (x)  of  the  initial  state,  is 
equivalent  to  its  convergence  to  some  mean  value  of  f(x )  as  x  is  averaged  over  all  possible 


19 


paths.  This  in  a  sense  provides  an  intuitive  justification  for  the  convergence  of  a  heat 
equation  to  a  constant.  The  case  in  point  arises  when  we  are  faced  with  a  deterministic 
diffusion  Xt,  which  is  alternatively  expressed  as  dXt  =  [10]dt.  The  corresponding  solution 
obtained  from  Eq.  (3.3)  is  U(x)  =  f(x0 1  +t,x0-2),  where  x0  is  the  initial  state,  clearly 
non-constant  as  expected. 

In  light  of  the  foregoing  development,  and  for  better  insight  and  intuitive  clarity,  we  find 
it  useful  to  carry  out  most  of  the  exposition  and  the  analysis  in  a  discrete  setting.  Prior 
to  delving  into  our  formulation  and  interpretation  of  a  Non-Linear  (NL)  diffusion,  we 
present  an  illustrative  example  where  a  so-called  controlled  diffusion  leads  to  a  Markov 
Chain  as  a  result  of  an  optimization. 


3.3.2  Finite  Markov  Chain  Example 

We  start  with  a  finite  Markov  chain  [31]  as  in  Fig.  3.2  with  three  possible  states  a,/?, 7 
with  respective  costs  1,  2,  3.  Our  goal  is  to  select  a  “best”  strategy  (minimum  cost)  for  a 
particle  to  make  a  two-step  transition  from  a  state,  say  a.  The  corresponding  transition 
probabilities  are  obtained  from  the  following  sets  of  strategies: 

Ga  =  {(l,0,|),(|.i.O)} 

G,  =  {(|.0,|),(|,i.O)} 

G-,  =  {(1,0,0),  (-.0, -),  (-.0, -)}. 

For  state  a  for  instance,  two  choices  are  possible: 

•  the  first  strategy  has  it  move  to  state  7  with  probability  and  remain  stationary  with 

probability  |, 

•  the  second  lets  it  move  to  state  /?  with  probability  \  and  remain  stationary  with  prob¬ 

ability  |. 

Taking  into  account  the  respective  costs  as  well  as  the  transition  probabilities,  an  optimal 
strategy  for  a  is  determined  to  be  (|,|,0)  with  a  minimum  cost  of  |,  while  in  state 
/ 3  a  cost  of  |  with  the  strategy  (|,|,0),  and  in  state  7  we  obtain  a  cost  of  1  with 
strategy  (1,0,0).  A  second  step  transition  may  be  similarly  found,  with  respective  costs 
for  a, /3,7,  of  |,|,|  resulting  from  (f,0,  |),  (f,  0,  |),  (f,  0,  §)  respectively.  This  process 
may  be  continued  indefinitely. 

Note  that  more  complex  strategies  are  possible  and  may  be  constructed  for  the  intermedi¬ 
ate  steps,  e.g.,  variable  strategies  along  the  steps,  etc..  When,  on  the  other  hand,  a  given 
strategy  set  only  depends  on  the  previous  state,  it  is  referred  to  as  a  Markov  strategy.  It 
is  also  clear  from  the  foregoing  example  that  the  resulting  process,  by  way  of  its  transition 
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In  an  attempt  to  solve  Eq.(  3.8)  and  using  its  stochastic  interpretation,  we  proceed  to 
write 


Ut(x)  =  EttS(f(XT)),  (3.10) 

where  /(•)  =  1/(0,  x)  is  the  initial  data.  In  particular,  we  obtain  one  step  transition  as 

Ut(x)  =  .  Et^{Et-T{f(XT))) 

=  J  Et-rM(Xr))p{T,x,y)dy 

=  J  Ut-r{y)PT(y\x)dy.  (3.11) 

Note  that  the  above  diffusion  is  a  backward  diffusion  which,  as  mentioned  in  Section  2,  is 
treated  as  a  forward  diffusion  for  the  sake  of  clarity.  The  probability  PT(y\x)  should  then 
be  interpreted  as  P(&_T  =  y|£r  =  x)  to  emphasize  the  backward  evolution  in  time/scale. 

3.4.1  Discrete  Time/Scale  Evolution 

In  discretizing  x  and  t.  we  account  for  a  reverse  time  evolution  by  relabeling  time  “t  —  r” 
by  1  and  “t  -  nr"  by  n,  hence  making  a  backward  diffusion  equation  look  more  like  a 
forward  diffusion.  We  denote  by  Un(x)  the  value  of  the  solution  at  time  step  nr  and 
location/state  x.  Eq.(  3.11)  can  then  be  written  in  the  form, 

Un+i(x)  =  Yyn(y)Pn{y\x).  (3.12) 

v 

Since  the  solution  to  Eq.(  3.7)  is  a  Gaussian  transition  density  function,  it  characterizes 
the  evolution  of  a  particle  along  a  Brownian  trajectory  starting  at  xo  and  time  t.  Using 
the  fact  that  a  limiting  process  of  a  random  walk  is  a  Brownian  motion,  we  may  compute 
the  solution  to  Eq.  (3.12)  at  any  desired  discrete  time/scale.  At  the  first  time  step  r  and 
for  a  1-D  case,  we  can  write 

Vi{x)  =  \  f(x-S)+1-f(x  +  S),  (3.13) 

scenario,  we  have 

=  — /(a?i  —  5,x  2)  +  -/(x  1  +  8,x  2)  +  -f(xi,x2  —  5)  +  -/(x  i,x2  +  5), 

(3.14) 

both  of  which  are  the  result  of  an  averaging  process.  More  generally,  we  can  write  the 
solution  to  the  linear  heat  equation  as  a  discrete  expectation  for  respectively  1-D  and  2-D 
as 

Un»(x)  =  lun(x-S)  +  ^Un(x  +  6),  (3.15) 


while  in  a  2-D 
Ui{xux2) 


22 


(3.16) 

(3.17) 


Un+i(xi,x2) 


^Un(x  —  5,  y)  +  ^Un(xl  +  6,x2)  +  ^Un(x1,x2  -  8)  + 
-Un(xi,X2  +  8) 


Due  to  the  underlying  random  walker  moving  to  its  neighbor  with  1/2  probability  in  1-D 
(and  to  its  four  nearest  neighbors  with  probability  1/4  in  2-D),  it  is  clear  that  the  linear 
evolution  will  indiscriminately  smooth  away  sharp  features  along  with  the  noise. 


3.4.2  A  Stochastic  View  of  Perona-Malik  Diffusion 

The  linear  stochastic  differential  equation  led  to  a  linear  diffusion  by  way  of  a  Laplacian 
as  its  corresponding  infinitesimal  generator.  Using  this  development  as  an  inspiration, 
together  with  its  discrete  stochastic  formulation  and  interpretation,  we  proceed  in  an 
analogous  manner  to  rewrite  the  P-M  equation  to  be  interpreted  as  a  particle-based 
diffusion. 

Proposition  2.  Based  on  a  particle  system  interpretation,  P-M  equation  may  rewritten 
as 


Un+i(x)  =  pn{x,x  +  8)Un(x  +  8)  +pn(x,x  -  8)Un(x  -  8), 

[1  -  pn(x,x  +  8) +pn(x,x  -  8}Un(x).  (3.18) 


Proof:  The  proof  follows  immediately  from  the  discretization  of  Eq.(  3.9)  and  rewriting 
1/2 g  (|  Un(x  ±  8)  -  Un{x )  |)  =  pn(x,x  ±  8)  =  pn{^n+ i  =  x  ±  8  |  £„  =  x)  to  denote  the 
transition  probability  of  a  Markov  chain  {£(.)}  to  move  from  state  x  to  state  x  ±  8. 

A  similar  expression  for  a  2-D  signal  (image)  may  be  written  as 
Un+1(x1,x2) 

=  pns^1(xl,x2)Un(x1  +  8 ,  x2)  +  pf^lUn{x\  -  5,  x2) 

+P1E~l(xi,x2)Un(xi,x2  +  8)  +p?y1Un(x1,x2  -  8) 

+  [1  -p%+1(xux2)  -p^+1(x1;x2)  ~pl+1{xux2) -p^1(xl,x2)]Un(x1,x2), 

(3.19) 


where 


Ps+1(xux2) 

Pn+\x  UX2) 
Pe+1(xux2) 

Pvvl{xux2) 


=  ^g(\Un(x1  +  8,x2) -Un(xi,x2)  \) , 

=  7^(1  Un(x  -  8,y)  -  Un(xuX2)  |) , 

4 

=  -g{\  Un(xux2  +  5)  -  Un{xi,x2)  |) , 
=  ^g(\Un{xux2- 8) -Un{xl,x2)  \) 
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represent  the  transition  probabilities  of  the  underlying  Markov  chain  £n,  i.e., 

Ps  (£n+i  =  (xi  +  5,  x2)  |  £n  =  (xi,x2))  =  Ps+1  fu )  ^2)  •  These  equations  are  intuitively  ap¬ 
pealing,  in  that  the  random  walk  of  a  particle/pixel  (or  the  diffusion)  takes  place  accord¬ 
ing  to  the  prevailing  one  sided  gradient  at  position  {x\,  x2)  in  any  of  the  four  directions.  At 
time  step  n  + 1,  a  south  (resp.  north,  east,  west)  moving  walk  takes  place  with  probabil¬ 
ity  P5-+1(a;i,  £2)  (resp.  p1}/1(xi,x2),  p1^1  (x\,  x2) ,  p^y1  (xi,  x2)) ,  and  the  particle  remains  in 
place  with  probability  p”+1  (a)  =  1  -pns+l  (aq ,  x2) -p^+1  (aq ,  x2)  ~Pe+1  (aq ,  x2 )  -p^t1  (aq  ,x2). 
It  is  clear  here  that  the  transition  probability  of  a  random  walk  is  determined  by  the  gra¬ 
dient,  thus,  the  resulting  diffusion  is  well  controlled.  This  is  in  sharp  contrast  to  the 
linear  diffusion  where  the  random  walk  invariably  takes  place  with  a  constant  probability 
of  1/4.  Note  that  while  the  derivation  of  an  exact  SDE  corresponding  to  P-M  equation 
as  an  infinitesimal  generator,  is  interesting  in  and  of  itself,  it  requires  a  more  complex 
system  of  particles  which  is  of  no  additional  insight  and  of  little  relevance  to  our  stated 
goal  in  this  paper. 

3.5  Two  Sided  Gradient-Driven  Diffusion 
3.5.1  A  variation  on  a  theme 

As  discussed  in  Section  2,  at  each  scale  of  our  analysis,  the  mean  value  of  the  process 
£/(•,•)  is  evaluated  as  a  result  of  a  non-homogeneous  random  walk  with  the  transition 
probability  controlled  by  the  underlying  process  at  the  previous  scale.  In  addition,  and  to 
avoid  potential  stability  problems,  we  ensure  that  the  probability  of  a  jump  of  a  particle 
(pixel)  farther  than  an  immediate  neighbor  is  zero,  which  effectively  emulates  a  continuous 
diffusion.  Furthermore,  we  ensure  that  there  always  be  a  one  step  transition  of  a  particle 
to  its  neighbors  to  avoid  a  slowdown  in  convergence  due  to  likely  stationary  states  [28]. 
We  adopt  this  paradigm  to  construct  a  non-homogeneous  Markov  chain  whose  transition 
probabilities  are  based  on  the  current  particle  states  and  their  functional  value.  This 
results  in  a  set  of  consecutive  transition  steps  through  scales,  each  in  a  sense,  defining  a 
new  random  process  with  a  new  probability  transition. 

While  the  goal  in  signal/image  processing  is  to  maximally  smooth  out  the  noise,  we  are  also 
keen  on  achieving  a  solution  that  is  as  faithful  as  possible  to  the  initial  underlying  signal. 
To  thus  help  better  localize  the  homogeneous  regions  together  with  their  boundaries,  we 
use  in  our  transition  dynamics  a  bidirectional  gradient-based  “probability  measure” .  (sub¬ 
gradient  in  continuous  space).  Using  the  Szokefalvi-Nagy’s  inequality [3 2],  to  optimize 
the  gradient  energy  (to  delineate  regions),  we  have  to  minimize  the  following  energy 
expression, 
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£((/„+,)  =  £f(EWi)) 

X 

=  J][([/„(a:  +  5)-[/„(^))(/7n+1(x)-C/n(a;-(5))]2 

X 

+  [(t/n(o:-5)-[/n(a;))(f/n+1(x)-C/„(x  +  5))]2 

where  t/n+i(x),  assumed  to  result  from  Eq.(  3.11),  is  written  as 

Un+i(x)  =  Pn+ i(x  -  8\x)Un{x  -8)  +  Pn+i(x  +  5\x)Un(x  +  8), 


(3.20) 


(3.21) 


with  Pn+1(x  —  <5|x)  +  Pn+ i(x  +  5|x)  =  1.  Minimizing  Eq.(  3.20)  entails  an  appropriate 
choice  of  a  probability  measure  as  follows, 

Theorem  1.  The  transition  probability  solving  Eq.(  3.20)  isvjiven  by 
Pn+i{x  -  5|x)  =  P{£n+ i  =  x-  <5|£n  =  x}  with 


Pn+  l(x  -  tf|x)  = 


\  Un{x +  8)  -Un{x)  |2 

Un(x  -8)-  Unix)  |2  +  |  Unix  +  8)-  Un{x)  |5 


(3.22) 


where  f/n+i(x)  satisfies  Eq.(  3.21). 


(See  Appendix  A  for  Proof).  ■ 

For  a  2-D  image,  we  denote  the  transition  probability  by  pns+lix i,  X2)  =  P{£n+ 1  =  (^1  + 
8,x2)\£,n  =  (xi,  X2)}  (similarly  for  other  probabilities)  and  obtain  the  following  expression 
for  the  transition  probability 


Ps+1(x  i,x2)  = 


<S  +  N  +  S  +  W' 


(3.23) 


where 


M  =  \Unix-8,y)-Unixux2)  I2 
<S  =  I  Unix  1  +  8,  X2)  -  Unix  1,  X2)  |2 
£  =  I  C4(xi, x2  +  8)  -  Unix l,X2)  I2 
w  =  \Un(x1,x2-8)-Unix1,x2)\2 

Using  the  above  transition  probability,  our  newly  proposed  diffusion  is  written  as 

Un+ i(®i,x2)  =  Unixi  4-  x2)p5+1(xi,x2)  +  Unix  —  5,  r/)p^+1(xi,x2) 

Unix  1,  X2  +  8)pnE+\x  1,  X2)  +  Unix,  ij  -  1  (®x ,  x2) .  (3.24) 
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challenges  to  simple  gradient-based  segmentation  and/or  linear  filtering,  and  this  is  for 
the  most  part  due  to  their  impulsive  nature.  The  results  of  such  approaches  usually  re¬ 
sults  in  visually  unpleasant  and  quantitatively  inaccurate  results  if  at  all.  In  Fig.  3.6,  the 
potential  for  complete  diffusion  rather  quickly  (for  an  inadequate  choice  of  the  threshold 
parameter  in  the  P-M  equation)  is  exhibited  whereas  and  as  demonstrated  above,  the  new 
approach  will  stabilize  with  no  parameter  adjustment.  For  a  well  known  stopping  time 
and  well  chosen  parameter,  P-M  approach  performs  quite  well,  this  may,  however,  may 
turn  out  to  be  a  limitation  for  a  number  of  real  applications. 

In  Fig.  3.7  an  enhancement/De-blurring-like  effect  using  our  algorithm  is  also  demon¬ 
strated  for  a  checker  board.  The  price  for  avoiding  the  explicit  knowledge  of  the  stopping 
time  using  our  approach,  is  the  arising  of  block  effects  which  although  common  to  many 
existing  techniques  is  a  drawback.  Although  this  may  be  fine  tuned  away  for  pure  De- 
noising  purposes,  we  are  currently  looking  into  techniques  specifically  addressing  such  an 
issue.  In  Figure  3.7,  a  De-blurring  example  is  shown,  demonstrating  the  capacity  of  the 
algorithm  to  enhance  edges  and  again  stabilize  at  staircase  frictions. 

For  establishing  a  more  quantitative  measure  of  performance  we  use  the  figures  in  Fig¬ 
ure  3.8.  A  pixel  deviation  is  computed  and  a  Monte  Carlo  simulation  is  conducted  to 
construct  an  error  rate  curve  which  is  consistent  with  our  visual  assessment  and  displayed 
in  Fig.  3.9. 
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3.7  Conclusion 


The  proposed  stochastic  interpretation  together  with  its  link  to  controlled  diffusion  are 
shown  to  not  only  explicate  existing  techniques  and  their  limitations,  but  to  also  provide 
sufficient  insight  to  develop  other  novel  physically  and  geometrically  driven  methodologies. 
We  have  also  succeeded  in  resolving  in  part  a  well  known  and  long  standing  problem  of 
unknown  stopping  criterion. 


Chapter  4 

Geometric  Stochastic  Flows 


4.1  Introduction  T 

Closely  related  techniques  to  the  above  anistropic  diffusions  are  those  of  curve  evolutions. 
As  we  elaborate  further  in  the  next  section,  the  so-called  geometric  heat  flow  may  be 
viewed  as  a  descendent  of  the  linear  heat  flow,  in  spite  of  the  fact  that  the  general  class 
of  curve  evolutions  itself  need  not  be  so.  Curve  evolution  has,  in  recent  years,  emerged 
as  an  important  application  of  partial  differential  equations  (PDE’s)  in  image  processing, 
computer  vision,  and  computer  graphics.  Their  applications  to  individual  curves  such  as 
in  edge-detection,  skeletonization,  and  shape  analysis,  but  have  also  been  extended  to 
their  simultaneous  action  on  the  level  sets  of  an  image  in  a  number  of  geometrically  based 
anisotropic  smoothing  algorithms.  A  detailed  account  of  these  methods  may  be  found  in 
[33,  34]  and  also  in,  e.g.  [7]. 

If  run  for  too  long,  however,  even  large  scale  features  will  be  destroyed.  The  reason  stems 
from  the  fact  that  as  the  geometric  heat  flow  shrinks  any  closed  curve,  the  curve  becomes 
more  and  more  circular  (elliptical  in  the  case  of  the  affine  flow)  and  will  eventually  collapse 
into  a  single  point  [35].  It  is  therefore  not  always  possible  to  preserve  desired  features  in 
the  shapes  of  objects  (corners  for  example)  if  too  much  evolution  is  required  to  remove 
a  significant  level  of  noise.  Furthermore,  it  is  not  well  understood  how  these  curvature- 
based  filters  are  affected  by  different  noise  distributions  and  when  this  sort  of  problem 
may  occur.  Much  like  the  anisotropic  diffusion  problem,  the  curve  evolution  technique  has 
largely  been  addressed  in  the  literature,  outside  [2,  36],  in  a  purely  deterministic  setting. 
In  this  part  of  our  research  effort,  we  provide  a  stochastic  formulation  of  the  geometric 
heat  equation  and  use  the  resulting  insights  to  develop  a  new  class  of  curvature-based 
flows  and  anisotropic  diffusion  filters  which  preserve  desired  features  in  the  shape  of  an 
object.  Under  these  new  flows,  evolving  curves  take  the  limiting  form  of  a  polygon  (see 
[37]  for  evolutions  of  polygons  related  to  the  geometric  and  affine  geometric  heat  flows). 
The  resulting  diffusion  models  may  therefore  be  applied  for  much  longer  periods  of  time 
without  distorting  the  shapes  of  polygonal  objects  in  the  image,  thereby  mitigating  the 
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We  may  write  these  directions  in  terms  of  the  first  derivatives  of  the  image  as 


_  (lij;,  tiy)  £  _  (  Uyt^x) 

1  \/ux2  +%2’  1/ ux2  +  Wy2  ’ 

Since  these  constitute  orthogonal  directions,  we  may  exploit  the  rotational  invariance 
of  the  Laplacian  operator  and  re-write  the  linear  heat  equation  in  terms  of  these  two 
variables: 

ut  =  V  •  (Vu)  =  +  um 

where  um  and  denote  second-order  directional  derivatives  in  the  directions  of  77  and  £ 
respectively.  It  is  possible  to  derive  the  following  expressions 


'U'x'Uxx  2UXUyUXy  H"  'U'y'LLyy 


U'TJT)  ” 

Ux+Ul 

= 

*  2iXlx'liy'lixy 

+  uluyy 1 

ul  +  ul 

(4.2) 

(4.3) 


By  subtracting  the  normal  diffusion  component  (4.2)  from  the  linear  heat  equation,  which 
diffuses  isotropically,  we  obtain  the  following  anisotropic  model,  which  diffuses  along  the 
boundaries  of  image  features  but  not  across  them 


ut  =  Ug  = 


UiyUXX  2uxuyuxy  -f"  u~.u 


yy 


K  +  Uy 


(4.4) 


We  may  obtain  this  same  equation  in  a  completely  different  and  much  more  geometric 
manner  by  specifying  the  evolution  of  each  level  curve  in  the  image.  Let  C  denote  a 
particular  iso-intensity  contour  which  we  will  deform  over  time  via  the  following  flow, 

Ct  =  Css  =  kN  (4.5) 

where  s  denotes  the  arclength  parameter,  k  the  Euclidean  curvature,  and  N  the  inward 
unit  normal.  Equation  (4.5),  referred  to  as  the  Geometric  Heat  Equation  (GHE),  is  well 
known  for  its  smoothing  properties.  It  has  been  shown  by  Grayson  [35]  that  any  closed, 
embedded  curve  evolving  according  to  (4.5)  will  convexify  and  smoothly  shrink  to  a  single 
point  in  finite  time,  becoming  more  and  more  circular  along  the  way.  This  flow  is  also 
referred  to  as  the  curve  shortening  flow  since  it  corresponds  to  the  gradient  (descent) 
evolution  of  the  arclength  functional.  See  [39,  40,  41]  for  a  more  extensive  discussion  of 
the  many  properties  associated  with  this  flow.  Because  the  evolution  speed  is  a  function 
of  the  curvature  at  each  point  on  a  curve,  this  flow  gives  rise  to  a  Euclidean  invariant 
scale  space  (see  [42,  43,  38])  in  which  finer  features  are  removed  first,  followed  by  coarser 
features,  as  the  curve  evolves.  A  related  flow,  based  upon  the  affine  geometry  of  the  curve, 
is  given  by  Ct  =  k1^N  and  shares  many  of  the  same  properties  as  the  curve  shortening 
flow  but  gives  rise  to  a  more  general  affine  invariant  scale  space  (see  [42,  44,  45]). 
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(«*(*,  X ),  Uy(t,  x  ))l\[ujt,  x  f  +  Uy{t,  X  )2.  It  follows,  9(ux(t,  x ),  uy(t,  x ))  =  tan 
Using  these  equations,  and  defining  an  operator  AGHE  of  the  form 

AGnEu(t,x)  -  sin 2  6(ux(t,x),uy(t,x))  uxx(t,x) 

—  2sin9(ux(t,x),uy(t,x))  cos9(ux(t,x),uy(t,x))  uxy(t,x ) 

+  cos2  9(ux(t,  x),uy(t,  x ))  uyy{t,x),  (4.7) 

the  geometric  heat  equation  (4.4)  can  be  re-written  as 

ut(t,x )  =  AGliEu(t,x ),  (4.8) 

u(0,x)  =  u0(x),  (4.9) 

where  u0(x)  is  the  initial  level  set  function.  We  next  show  how  an  evolution  equation  in 
fact  corresponds  to  an  infinitesimal  generator  of  a  Stochastic  Differential  Equation  (SDE), 
by  using  Ito  diffusions  and  the  Kolmogorov  backward  diffusi^p  theorem  [46,  2,  36]. 

4.3.1  Ito  Diffusion 

The  diffusion  of  a  particle  is  usually  well  modeled  by  an  SDE  which,  in  turn,  represents 
the  underlying  microscopic  process  of  an  evolution  of  a  pixel  or  a  point.  The  dynamics  of 
this  evolution  at  a  macroscopic  level  are  captured  by  a  PDE,  henceforth  also  referred  to 
as  a  generator  (infinitesimal)  of  the  diffusion. 

Definition  2.  Suppose  we  want  to  describe  the  motion  of  a  small  particle  suspended  in  a 
moving  liquid,  subject  to  random  molecular  bombardments.  If  b(t,x)  E  Rn  is  the  velocity 
of  the  fluid  at  a  point  x  E  Rn  and  time  t  E  R+ ,  then  a  widely  used  mathematical  model 
for  the  position  X  (t)  of  the  particle  at  time  t  is  an  SDE  of  the  form 

dX(t)=b(t,X(t))dt  +  cr(t,X{t))dB(t),  (4.10) 

where  X  (t)  is  an  n- dimensional  stochastic  process,  cr(t,x )  G  R'lXm,  and  B  (t)  is  an 
m- dimensional  Brownian  motion,  b  (•,  •)  is  called  the  drift  coefficient,  and  cr  (-,  •)  is  called 
the  diffusion  coefficient. 

The  first  term  in  this  equation  corresponds  to  a  non-random/deterministic  motion,  whereas 
the  second  term  models  randomness  or  noise  in  the  motion. 

The  solution  of  such  an  SDE  may  be  thought  of  as  a  mathematical  description  of  the 
motion  of  a  small  particle  in  a  moving  fluid,  and  such  stochastic  processes  are  called  (Ito) 
diffusions  [46].  For  many  applications,  a  second  order  partial  differential  operator  A  can 
be  associated  to  an  Ito  diffusion  X  (t)  given  by  Eq.  (4.10).  The  basic  connection  between 
A  and  X  ( t )  is  that  A  is  the  generator  of  the  process  X  (t).  If  w(x )  E  Cl (Rra),  (i.e.,  it  is 
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continuous  with  continuous  derivatives  up  to  order  2,  and  has  a  compact  support),  then 
A  is  given  in  the  form 


Aw=  2 


J2(cr<TT),,(x) 


d2w 
dxi  dxj 


+  ^6i(a:) 
i 


dw 

dxi' 


(4.11) 


Next,  we  state  a  theorem,  the  so-called  Kolmogorov’s  backward  equation  [46],  which  gives 
a  probabilistic  solution  to  linear  PDE’s. 

Theorem  2.  Definey(t,x)  =  Ex[f(X  (t))],  where  X  (t)  =  X^(t)),  andEx[-] 

is  the  expectation  operator  with  respect  to  the  probability  law  of  X  ( t )  starting  at  the  point 
x ,  then  there  exists  an  operator  A  such  that, 


dy 

dt 

7(0,*) 


=  Ay,  t  >  0,  x  Gl2, 
=  f(x ),  x  G  R2. 


(4.12) 

(4.13) 


4.3.2  Stochastic  Formulation 


T 


In  light  of  the  foregoing  development,  a  natural  question  which  arises  is:  given  a  PDE 
which  governs  a  curve  shortening  flow,  can  we  obtain  a  corresponding  SDE  associated 
with  the  underlying  diffusion? 

The  nonlinearity  of  GHE  presents  a  significant  challenge  to  find  a  global  Ito  diffusion  which 
explains  the  overall  microscopical  behavior  of  the  system.  Our  approach  here  for  solving 
such  a  nonlinear  problem  is,  to  explore  the  short-time  behavior  by  linearizing  around  a 
known  (nominal)  solution.  The  perturbation  equations  so  obtained  will  be  linear  and 
hence  an  approximate  solution  to  the  nonlinear  problem  can  be  obtained  as  the  nominal 
value  plus  the  perturbation  term.  Let  us  denote  by  un(t,x)  the  solution  to  Eq  (4.8): 

=  sin2  (tan-1(— )^1  —  sin  f2tan-1(— )^  uf.  +  cos2  ftan-1(-^) 

dt  V  X7  V  X  /  y  V  X 

and  if  we  write  u(t,x)  as 

u(t,  x )  =  un(t,  x )  4-  eu(t,  x ), 

un(t  X ) 

and  define  the  corresponding  nominal  angle  9n(t,  x )  =  tan~1(tt^f’a.^),  we  get  a  linearized 
version  of  the  geometric  heat  equation  around  a  nominal  value: 

du^X).  «  ^lGHElinu(t,  x )  =  sin2(0n(x ))  uxx(t,  x )  -  sin(20,l(£c ))  uxy(t,  x )  + 

cos2(6*n(a; ))  uyy(t,  x )  +  c(a? )(-«”( x)  ux(t,x )  +  v%(x)  uy(t,x)), 

(4.14) 

where  c(x )  =  [sin(20"(a:  ))«x(a! )  ~  «»(* ))  -  cos{2dn(x  ))2 unxy{x )]  (see 

Appendix  B  A  for  details  of  this  derivation). 

In  light  of  this,  we  can  proceed  to  state  the  following: 
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Proposition  3.  The  right  hand  side  of  the  linear  PDE  in  Eq.  (f.lf)  is  the  generator  of 
the  following  Ito  diffusion  satisfying  the  SEE 


(  dX  W(£) 

V  dXW\t) 


=  c(X(t)) 


~uny(X  (t)) 
<(X  (t)) 


dtp  y/2 


-sin 6n{X  (t)) 
cos  6n{X  (t)) 


dB(f), 

(4.15) 


Proof.  The  operator  *4GHEiin  in  Eq.  (4.14)  is  first  re-written  as, 


Ashbu.  =  br(I)-V+-(T(X)<TT(I)0H 


=  c(X) 


-k(x)Y.(£ 


K(x) 


t  + 

dy  ) 


sin2  0n(X )  —  sin  6n  (X )  cos  0n {X 

—  sin  9n  (X  )  cos  6n  (X )  cos20n(X) 


©H, 


where  H  is  a  Hessian  operator  and  0  is  a  Hadamard  product.  The  factorization  of  \cr  crT 
leads  to 


<r(X)  =  V 2 


-sinr(x; 
cos  9n(X ) 


(  -un(X 

and  by  identification,  b  (X )  =  c(X )  nyv . 

\  ) 


Given  the  functions  b  (X  ( t )),  and  <r  (X  (f)),  we  come  up  with  a  pair  of  processes  (X  (t),  B(t)) 
such  that  the  SDE  in  Eq.  (4.15)  holds.  In  this  case,  the  solution  X  ( t )  is  called  a  weak 
solution,  as  it  does  not  specify  beforehand  the  explicit  representation  of  the  white  noise, 
i.e.  the  version  B(t)  of  the  Brownian  motion  is  not  given  in  advance.  □ 


Both  the  drift  and  diffusion  coefficient  vectors  of  this  SDE  are  in  the  tangent  direction 
of  our  level  curves,  which  helps  us  interpret  it  as  a  1-dimensional  Ito  diffusion  on  the 
instantaneous  tangent  direction  T  (u2(t),Uy(t)).  A  differentiability  assumption  on  u(t,  x ) 


u(t  +  5t,  x )  —  u(t,  x )  duit,x) 

hrn - — - - - -  =  — 5- — 

5t — >-0  oi  dt 


AotmiMt*  x ) 


is  sufficient  for  a  short-time  existence  of  the  linearized  PDE  version  of  the  nonlinear 
geometric  heat  equation. 

Using  Kolmogorov’s  theorem  cited  above,  and  assuming  u(t,  x )  and  its  derivatives  are 
“sufficiently  regular”  (Lipschitz  properties),  starting  at  each  time  t,  the  diffusion  X  ( t ) 
in  Eq.  (4.15)  is  constructed  for  each  time  interval  ( t  —  8t,t ),  and  may  be  used  to  write  a 
Backward  Kolmogorov  Equation, 


u(t  —  8t,x)  =  E{u(t,X  (i))/X  ( t  —  St)  =  x  }, 


as  a  mean  value  around  each  pixel  dictated  by  the  motion  of  the  constructed  diffusion 
X  ( t ).  This  equation  can  also  be  written  in  forward  time  (since  in  the  small  time  step 
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not  change  but  only  the  measure  on  the  trajectories  change.  This  theorem  involving  a 
change  of  measure  provides  us  with  a  means  of  changing  the  mean  of  the  process  X  ( t ) 
we  obtained  in  Eq.  (4.16),  particularly  removing  the  drift  and  obtaining  the  process  in 
Eq.  (4.17). 

This  intuitively  appealing  interpretation  of  a  particle/pixel  motion  in  the  process  of  a 
diffusion  is  shown  in  the  next  section  to  be  particularly  useful  and  insightful  for  developing 
more  general  and  feature/shape  adapted  flows. 


4.4  A  New  Class  of  Flows 

The  insight  gained  from  the  .tangential  Brownian  motion  on  a  curve  together  with  the 
normal  angle  8(t,x),  lead  to  the  idea  of  constraining  the  Brownian  motion  at  some 
specific  orientation  angles  at  each  point  x .  A  natural  generalization  of  the  geometric  heat 
equation,  based  upon  the  stochastic  framework  presented  infection  4.3,  is  to  construct 
an  SDE  weighted  by  a  carefully  chosen  functional  h(8n),  (h(.)  €  C°°(Rn))  designed  to 
capture  specific  features  in  an  image,  and  we  write  locally 

dX  (t)  =  Tn(X  (t))  {cn(X  (: t))h(8n{X  (*)))  dt  +  V 2  h(8n{X  (t)))  dB{t)). 

Here,  neglecting  the  drift  motion  concentrating  on  pure  diffusion,  the  Brownian  motion 
in  the  tangent  direction  is  being  further  constrained  at  some  specific  orientation  values, 
i.e.  at  the  zeros  of  the  h(6n)  function, 


dX  ( t )  »  y/2  Tn(X  (t))  h(8n(X  (t)))  dB(t)).  (4.1.9) 

Constraining  the  diffusion  of  particles  at  points  with  specified  orientations  is  aimed  at 
extracting  desired  features  of  a  contour  as  it  is  being  smoothed.  Such  models  are  generated 
by  the  following  class  of  PDE’s  which  generalize  the  geometric  heat  flow  (4.4) 

^LL  =  h?(e(t,  x))um,x),  (4.20) 

which  is  locally  the  generator  of  the  diffusion  in  SDE  (4.19).  When  applied  to  an  image, 
this  flow  induces  the  following  curve  evolution  on  each  iso-intensity  contour  C 

r)C 

—  =  ti2(6)KN.  (4.21) 

4.4.1  Well-posedness  of  the  generalized  model 

Proposition  4.  The  corresponding  PDE’s  (4-20)  are  well-posed. 
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Proof.  The  geometric  heat  equation  which  corresponds  to  the  simplest  case  of  this  class 
with  h2(9(t,x))  =  1,  Vf,  Mx.  has  been  shown  to  be  well-posed,  and  its  existence  and 
uniqueness  properties  may  be  found  in  [43,  47,  48]. 

The  operator  of  the  geometric  heat  equation  is  given  by 

L[«]  =  £[«]-  ^  =  0  (4.22) 

where 

sin2  9  uxx  —  2  sin  9  cos  9  uxy  +  cos2  9  uyy,  (4.23) 

is  the  principal  part  of  the  operator  L.  The  matrix  of  coefficients  [a;j]  is  positive  semi- 
definite  with  the  eigen  values  1  and  0.  If  we  multiply  this  matrix  by  a  positive  function, 
it  remains  positive  semi-definite.  Such  elliptic-parabolic  operators  satisfy  a  maximum 
principle  (see,  for  example,  [49]).  In  our  case  we  multiply "%  a  non-negative  function 
h2(9)  which  can  be  made  strictly  positive  by  adding  a  very  small  number,  e  >  0, 

[h2(9)  +  e\Cu  >0. 

This  results  in  a  family  of  nonlinear  parabolic  equations  each  of  which  satisfies  a  strong 
maximum  principle.  Our  operator  is  obtained  in  the  limit  as  e  — »  0.  □ 


4.4.2  Polygon  leading  diffusions 


The  geometric  heat  equation  is  a  rotationally  invariant  flow  which  evolves,  as  mentioned 
earlier,  any  shape  into  a  circle  [35].  It  is  the  only  rotationally  invariant  shape  in  Euclidean 
space.  If  we  wish  to  capture  more  general  shapes  (triangles,  .squares,  etc  ...  )  it  is  only 
then  natural  to  consider  flows  which  are  not  rotationally  invariant.  Such  a  class  is  given 
by  the  form  (4.21)  when  h{9)  is  chosen  to  be  other  than  a  constant.  If  we  are  particularly 
interested  in  polygons,  we  may  consider  periodic  functions  (whose  periodicity  is  dictated 
by  the  desired  number  of  vertices)  such  as 


h2{9) 


cos2  ( n9 ) 
sin2  ( n9 ), 


(4.24) 


leading  to  curve  evolution  equations  of  the  form 

f)C  f)r 

—  =  cos2 {u9)kN  or  —  =  sin2(n0)/clV  .  (4-25) 

If  we  apply  (4.25)  to  a  convex  shape,  there  will  be  In  points  on  the  curve  which  do  not 
diffuse  (corresponding  to  the  zeros  of  cos(n0)  or  sin (n9))  at  equally  separated  rotations 
of  the  unit  normal  N .  As  the  unit  normal  moves  further  and  further  away  from  these 
angles,  the  diffusion  increases.  It  hence  makes  sense  to  expect  a  curve  to  develop  vertices 
(points  of  maximal  curvature)  at  these  points. 
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Lemma  1.  The  angle  of  a  unit  normal  does  not  change  at  points  where  the  chosen  func¬ 
tion  h2(6)  vanishes.  Those  points,  in  turn,  remain  fixed  for  a  short-time,  and  their  speed 
remains  at  zero. 


Proof.  Assume  that  a  family  of  curves  C(t,p),  where  p  is  any  parameter  along  the  curve, 
evolves  according  to  the  evolution  equation 

dC 

—  =  a(t,p)T  +p(t,p)N  (4.26) 

The  evolution  equation  for  the  angle  of  the  unit  normal  is  given  in  [39]  as 

99  -1  ta  1 

where  g  —  ||CP||  =  sjXp  +  Tp  is  the  length  along  the  curve  (metric).  If  we  consider  the 
case  a  =  0  and  (3  =  —h2{8)n  (following  the  convention  used  by  the  authors  in  [39]), 
which  corresponds  to  the  form  of  the  deformation  we  proposed,  the  orientation  evolution 
is  governed  by 

I  =  \  ^ 

=  i  {2  h{8)  (h(6))p  K  +  h2(9)  kp) 

=  i  h{8)  {2  {h(8))p  k  +  h{8)  kp}  ,  (4.27) 

Notice  that  %  —  0  for  those  points  at  which  h(8)  =  0.  □ 


We  note  that  in  [39],  the  orientation  of  a  curve  is  defined  as  the  angle  subtended  by  the 
tangent  and  the  x-axis,  whereas  we  here  define  8  as  the  angle  subtended  by  the  normal 
and  the  x-axis.  There  is,  however,  complete  equivalence  in  so  far  as  the  evolution  equation 
of  the  angle  8  is  concerned. 


In  light  of  the  above  development,  we  can  thus  state  that  the  zeros  of  the  function  h{8) 
lead  to  fixed  end  points  of  curve  segments.  Fixing  two  end  points,  say  Oi  and  a2  ,  we 
examine  the  evolution  of  curvature,  whose  general  form  is  given  by  (in  [39]) 


3  k  d2P 

dt  ds 2 


5 


where  s  is  the  arc-length  parameter  along  the  curve.  When  substituting  a  =  0  and 
/ 3  =  —h2(8)t c  into  this  equation,  we  have 


ft  ■=  [h2(8)K]ss  +  h2(8)K3 

r)  K 

—  =  [(h2(8))ssK  +  2(h2(8))sKs  +  1i2(8)kss]  -I-  h2(8)n3 

—  fi2(9)  Kgs,  +  hf{8)  ac3  +  (h2(fl))ss  «  +  2(h2(0))s 

diffusion  term  reaction  term 
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This  clearly  demonstrates  that  a  regularizing  diffusion  takes  place,  since  the  multiplica¬ 
tive  factor  h2(9)  never  becomes  negative  (which  would  result  in  an  ill-posed  backward 
diffusion).  In  addition,  we  have  the  reaction  term  which  is  composed  of  functions  of  k, 
/v3,  and  ks. 

We  have  hence  shown  that  with  fixed  end  points,  a  particular  curve  segment  subject  to 
the  new  evolution  equation  for  the  curvature  shown  above,  results  in  a  straight  line  as  a 
final  solution. 

Now,  we  can  state  a  theorem  where  we  put  our  argument  of  convergence  to  regular 
polygons. 

Theorem  3.  A  convex  curve  C  subject  to  the  evolution  Ct  =  Ii2(9)kN  will  converge  to 
an  M-sidedj  regular  polygon  whose  M  vertices  will  be  formed  at  those  vanishing  points  of 
the  function  h2(9). 


The  proof  of  this  theorem  can  be  completed  using  the  arc-length  evolution  equation 


where  ds  denotes  the  incremental  arclength  of  C.  Since  the  integrand  is  strictly  positive,  we 
see  that  a  curve  will  continue  to  shrink  until  curvature  vanishes,  that  is  the  curve  segment 
converges  to  a  straight  line  between  the  end  points  a\  and  a2.  This  in  conjunction  with 
the  above  lemma  completes  the  proof  of  the  theorem. 


4.5  Experimental  Results 

4.5.1  Examples  in  Polygonization 

To  substantiate  the  stated  theorem,  we  next  present  examples  illustrated  in  Fig.  B.5  and 
Fig.  B.6.  In  our  experiments  with  contours,  we  use  the  narrow-band  implementation  of  the 
level  set  method  developed  in  [50].  The  time  step  is,  8t  =  0.2.  Numerical  implementations 
of  the  proposed  flows  are  applied  to  a  variety  of  convex  shapes  shown  in  part  (a)  of  each 
figure.  In  Fig.  B.5,  the  shapes  in  part  (b)  were  obtained  by  using  h2(9)  —  cos2  (29)  via 
the  following  curve  evolution 


Ct  =  cos2  (2 9)kN  , .  (4.28) 

while  the  shapes  in  part  (c)  were  obtained  using  h2(9)  =  sin2  (29).  In  both  cases,  we 
expect  to  obtain  four-sided  and' regular  polygons.  The  zeros  of  cos(2$)  and  the  zeros 
of  sin (29)  are  however  45  degrees  out  of  phase.  As  such,  we  see  the  evolved  shapes 
in  part  (b)  taking  the  form  of  a  square,  whereas  the  evolved  shapes  in  part  (c)  take 
the  form  of  a  diamond,  corresponding  to  a  45  degree  rotation  of  the  shapes  in  part 
(b).  In  Fig.  B.6,  we  see  the  effect  of  using  different  periods.  The  shapes  in  part  (b) 
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are  obtained  using  Ct  —  cos2  (3(9) kN  ,  while  the  shapes  in  part  (c)  are  obtained  using 
Ct  =  sin2(1.5((9  —  tt/2))kN  .  In  the  first  case,  we  expect  6  vertices,  and  in  the  second 
case  we  expect  3  vertices.  Our  expectations  match  the  results  shown  in  part  (b)  and  (c), 
where  we  observe  hexagonal  and  triangular  shapes,  respectively. 

We  may  also  apply  these  flows  to  the  level  sets  of  an  image  in  the  same  manner  that 
the  geometric  heat  equation  may  be  applied.  This  gives  rise  to  a  family  of  anisotropic 
smoothing  filters  which,  unlike  the  geometric  heat  equation,  are  not  rotationally  invari¬ 
ant.  This  feature  can  be  useful  in  smoothing  noisy  images  where  corners  and  edges  are 
priorly  known  to  have  certain  orientations.  These  diffusions  are  modeled  by  PDE’s  of  the 
following  form: 

u,  =  h2(0)V  -  ||V«||.  (4.29) 

Note  that  the  trigonometric  expressions  we  have  considered  for  h2(9)  can  be  written  in 
terms  of  the  first  derivatives  of  u,  for  example 


cos2 (2 9)  = 


Ox  +  Ul)2 


and 


sin2  (29)  = 


(2  llx'lly) 

{ul  +  u2)2  ’ 


allowing  one  to  implement  the  PDE  without  having  to  compute  the  orientation  of  the  unit 
normal  to  each  level  curve.  Note  that  these  expressions  involve  only  first  order  derivatives 
and  therefore  do  not  alter  the  quasi-linear  structure  of  these  second  order  flows. 


4.5.2  Examples  in  Feature-Preservation 

Feature-preserving  properties  as  well  as  polygonal  approximation  properties  of  the  pro¬ 
posed  flows  will  be  demonstrated  in  this  section.  We  illustrate  the  idea  of  capturing 
different  polygonal  features  of  shapes  by  our  proposed  flows  on  the  following  examples. 

The  first  example  is  a  “chef”  shape  with  both  round  and  polygonal  features  as  shown 
in  Fig.  B.7.  Geometric  heat  flow  Ct  =  kN  ,  (n  =  0),  evolves  these  features  into  circles 
as  shown  in  the  second  row  of  Fig.  B.7  for  time  points  t  =  40,80,160.  Particularly,  at 
t  =  160,  most  parts  of  the  shape  turns  into  incomprehensible  blob-like  structures.  In 
contrast  to  this,  polygonal  features  of  the  chef  like  his  nose,  and  tray,  are  preserved  by 
the  flow  Ct  =  sin2(29)hiN ,  (n  =  2),  which  favors  diamond-like  structures  (see  third  row 
of  Fig.  B.7  for  t  =  40,80,160).  Similarly,  the  flow  Ct  =  cos2(A9)k,N  ,  (n  =  4),  favors 
octagonal  features  as  shown  on  the  fourth  row  of  Fig.  B.7,  which  is  observed  at  chef’s 
hat  at  all  time  points  t  =  40,80, 160.  The  regularity  of  these  flows  is  readily  observed 
through  the  smoothness  of  the  resulting  shapes.  When  we  view  each  row  from  left  to 
right,  we  observe  a  progression  from  finer  to  coarser  scale.  The  scale-spaces  produced 
by  our  modified  flows  in  the  last  two  rows  are  visually  more  pleasing  since  corners  are 
preserved,  whereas  in  the  row  above  we  see  them  smoothed  away  by  the  pure  geometric 
heat  flow. 
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The  second  shape  example  is  a  fish  which  contains  some  fine  detail  structures  as  well  as 
coarse  features  (Fig  B.8).  The  second  row  shows  the  result  of  the  geometric  heat  flow 
Ct  =  kN  ,(n  =  0),  which  smoothes  away  not  only  fine  features  but  some  coarse  features 
as  well  (the  fins  for  example).  The  results  of  the  flow,  Ct  =  cos2  (29)  kN  ,(71  =  2),  are 
shown  in  the  third  row  of  Fig  B.8.  In  this  case,  rectangular  features  are  preserved  for 
longer  periods  throughout  the  evolution.  Finally,  the  flow,  Ct  =  sin2(A9)KN ,  (n  =  4), 
is  depicted  in  the  last  row,  preserving  octagonal  features  as  shown  in  the  nose  and  the 
dorsal  things. 

In  the  third  example,  we  start  with  a  noisy  shape  at  time  t  —  0  shown  in  Fig.  B.9. 
This  shape  is  evolved  with  the  geometric  heat  flow  Ct  =  k N ,  (n  =  0) ,  the  flow  Ct  = 
sin2(1.5(0  -  7r/2 ))kN  ,  (n  =  1.5),  the  flow  Ct  =  sin2 (29) kN  ,  (n  =  2),  and  the  flow 
Ct  =  cos2(2.5(0  —  -k/2))kN  ,  (?r  =  2.5),  as  shown  in  Fig.  B.9.  The  geometric  heat  flow 
at  the  top  row  quickly  smoothes  corners  of  the  shape  out,  and  at  coarser  scales,  the 
shape  loses  all  of  its  features.  The  initial  shape  converges  to  a  circle  regardless  of  the 
global  feature  of  the  plane  being  a  polygonal  shape.  This  motivates  the  application  of  the 
geometric  heat  flow  with  a  sin2  (nO)  factor,  where  n=1.5  whole  weak  limiting  shape  is  a 
triangle  which  intuitively  matches  the  coarser  form  of  the  given  plane  shape.  Similarly, 
for  n  =  2  and  n  —  2.5,  different  features  of  the  shape  are  kept,  and  persist  over  a  much 
longer  time  period  as  can  be  observed  from  the  column  of  shapes  at  t  =  400.  Note  that 
the  geometric  heat  flow  result  at  the  top  quickly  washes  out  any  similarity  to  the  actual 
shape,  whereas  the  results  of  the  other  three  flows  keep  the  global  shape  as  well  as  some 
finer  details  on  the  wings,  the  tail,  and  the  head  part. 


4.5.3  Examples  with  Gray-scale  Images 

The  proposed  flows  may  also  be  applied  to  images  in  a  straightforward  fashion.  For  the 
case  h2(9)  =  cos2  (29),  all  level  sets  of  the  image  are  driven  to  rectangles,  thereby  enhanc¬ 
ing  those  features  in  an  image.  Such  features  can  be  found  in  contemporary  buildings 
where  one  example  in  NCSU,  Centennial  Campus  is  shown  in  Fig.B. 10(a).  The  part  of 
the  building  image  with  an  additive  Gaussian  noise  is  shown  in  Fig.B. 10(b).  The  2nd  row 
shows  the  results  of  the  geometric  heat  flow  ut  =  ug  at  t,  —  10, 20,  50.  The  noisy  image  at 
t  —  0  is  smoothed  out  very  quickly  at  the  expense  of  rounding  off  all  the  corners  because 
the  level  sets  of  the  image  converge  to  circles.  The  3rd  row  shows  the  ut  =  cos2(29)u # 
flow  results  at  the  same  time  points  t,  =  10,20,50.  Since  the  diffusion  is  constrained  in 
order  to  drive  image  level  sets  to  rectangles,  the  removal  of  noise  is  slower.  However, 
the  rectangular  features  still  nicely  appearing  after  noise  removal  (see  the  image  on  the 
right),  makes  it  worthwhile  to  slow'  down  the  De-noising  effect  of  the  geometric  heat  flow 
according  to  our  needs. 

In  Fig.  B.ll,  an  experiment  involving  diamond-like  shapes  in  the  image  taken  from  a 
poster  on  a  wall  is  shown.  In  the  middle  row,  rounding  effects  on  diamond  shapes  in 
this  image  performed  by  geometric  heat  flow  is  clearly  observed  during  evolution.  The 
proposed  flow,  shown  in  the  bottom  row,  which  takes  the  form  ut  =  sin2(20)u^  for  this 
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particular  shape  is  very  good  in  carrying  out  a  shape-adapted  smoothing  which  takes 
place  at  the  boundaries  of  the  diamonds.  The  slight  blurring  effect  on  the  picture  at 
continued  application  however  is  due  to  the  interaction  between  consecutive  level  curves. 

A  photograph  taken  by  path-finder  in  mars,  shown  in  Fig.  B.12,  is  argued  to  be  a  hexagon¬ 
shaped  structure  on  mars’  surface.  The  particular  flow  adapted  to  this  shape  is  given  by 
ut  —  cos2 (3 and  the  resulting  images  at  the  second  column  of  the  figure  demonstrate 
a  better  smoothing  performance  at  the  boundaries  of  the  hexagon  when  compared  to  the 
images  in  the  first  column  processed  by  geometric  heat  equation.  From  low  scales  to  very 
large  scales,  the  hexagon-adapted  flow  enhances  and  keeps  on  highlighting  the  related 
structure. 

A  noise  contaminated  Aerial  image  is  shown  in  Fig.  B.  13 (a).  The  geometric  heat  equation 
(see  2nd  row,  t  —  20, 40, 80)  sweeps  away  the  shape  information  of  the  important  details 
such  as  the  city  on  the  left  bottom,  the  white  bright  rectangle  on  the  right  bottom,  and 
the  black  feature  at  the  top.  The  three  images  resulting  from  the  Ut  =  cos2 (20)  flow, 
are  quite  sharp  at  the  edges  between  both  low  and  high  contrast  fields,  therefore  more 
useful  in  recognition  of  details  as  well  as  removing  noise. 

A  last  example  is  shown  in  Fig.  B.14,  where  windows  and  a  roof  of  a  section  of  a  house 
are  seen.  On  the  left  column,  are  the  results  of  the  geometric  heat  flow  ut  =  at  times 
t  =  40,80, 160,  and  on  the  right  are  the  results  of  ut  =  cos2(29)u^  at  same  time  steps. 
The  noise  is  successfully  removed  by  the  geometric  heat  equation  whose  smearing  effect 
on  different  regions  into  one  another  is  also  slow,  at  a  cost  of  a  problematic  rounding  off 
of  corners.  At  time  t  =  160  for  the  result  on  the  right  bottom,  approximately  the  same 
amount  of  noise  as  that  of  geometric  heat  equation  at  t  =  40  is  removed,  and  in  addition 
to  that  the  corners  are  still  well-preserved. 

4.6  Conclusions 

Using  insights  from  a  stochastic  formulation  of  the  geometric  heat  equation,  we  have  pre¬ 
sented  an  alternative  stochastic  view  of  the  nonlinear  filtering  with  an  ability  to  generalize 
and  propose  new  flows  driven  by  desired  geometrical  features,  such  as  polygonal  shapes 
of  interest  in  a  variety  of  shape  recognition  tasks. 
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Appendix  A 


Proof:  We  first  proceed  to  re-express  £(Un+i)  in  terms  of  Eq.(  3.21)  and  Pn+i(x  -1-  d|x)  = 
1  —  Pn+i(x  —  <5|.r).  By  subsequently  differentiating  with  respect  to  Pn(x  —  8\x)  and  bearing 
in  mind  that  a  two  sided-gradient  is  used,  we  have  the  following  equation 

d(£(Un+1))/dPn+1(x-8\x)  T 

=  2[(Un(x  +  8)-Un(x))2(Un+l(x)-Un(x-8))](dUn+i(x))/dPn+1(x-8\x) 

+  2 [(Un{x  ~8)~  Un(x))2(Un+\(x)  -  Un{x  +  8))}(dUn+1(x))/dPn+i(x  -  6\x) 

(A.l) 

Setting  d{£{Un+i))/dPn+i{x  —  <5|a*)  =  0  ,  and  assuming  the  non-degenerate  case  of 
dUn+i(x)/dPn+i(x  —  <5 1 a;)  ^  0,  the  optimal  transition  probability  at  the  n-th  step  im¬ 
plies 


[(Un{x  +  8)-  Un(x))2  +  (Un(x  ~  6)  ~  Un(x))2]Un+1(x) 

-  (Unix  +  5)-  Un(x))2Un(x  -8)  +  (Unix  -8)-  Un(x))2Un(x  +  <5)  (A. 2) 

where  the  replacement  of  Un+\  (x)  with  Eq.(  3.21)  will  reduce  to  Eq.  (  3.22).  Note  that  if 
dUn+\(x) / dPn+\(x  —  8\x)  —  0,  we  can  see  that  a  left  sided-gradient  is  equal  to  the  right 
sided-gradient  resulting  in  an  optimal  choice  of  probability  1/2.  ■ 
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Appendix  B 


Let  us  denote  by  un(t,x )  the  solution  to  Eq  (4.8): 


If  we  write  u(t,  x )  as 


then  eu(t,  x )  satisfies 


u(t,  x )  =  un{t,  x )  +  eu(t,  x ) 


deu 

dt 


=  «  +  eux,  Uy  +  euy,  unxx  +  euxx)  -  fi(ux,uy,  unxx)  -  {/2«  +  eux, uny  +  euy,uxy  +  euxy) 


~h(uxi  uy->  uxy)  I  +  h(ux  +  euxi  Uy  +  euy>  uyy  +  euyy)  h(ux>  uy>  uyy)' 


(B.l) 


Assuming  /i(-,  •,  •),  /2(*,  •,  •)  and  /3(-,  •,  •)  are  differentiable  in  their  arguments,  we  can  ex¬ 
pand  f\ (•,  •,  •)  in  Taylor  series  about  (ux,  uy,  uxx),  /2(-,  •,  •)  about  («",  uy,  uxy),  and  /3(-,  •,  •) 
about  ).  For  notational  simplicity,  let  us  denote  by  xx  -  (ux,uy,uxx)  and 

sc?  =  («;,u”,«;x),  s2  =  {ux,uy,uxy)  and  =  {unx,u*,unxy),  and  x3  =  («*,«»,  «w)  and 
x”  =  {ux,Uy,Uyy).  If  we  assume  that  eu(t,x)  is  small  enough,  we  can  neglect  higher 
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order  terms  and  write  a  linear  approximation  as 

^  ^  ex i  •  V/i(*“ )  -  ex2  ■  V/2(*J )  +  ex3  •  V/3(a:£ ) 


(B.2) 


/ 


<9eu 

~dt 


€('UX  Wxx) 


2 sin  (tan  K^f))  cos  (tan 
2  sin  (tan-1(g))  cos  (tan-1(g))  ^ 
^  sin2  (tan-1(g 


i+(^)2  («5)2  Uxx 

-^\un 


*  \ 


<>2  U?  XX 


c(^.T  ^xy) 


2  cos  (2  tan  1(^4)'j  — 

V  Ku*’)  i+c^ft 


_ qtn 

vn 

n\  2 


2  cos  (2  tan-1  (-£)')  — ^r-ipr 


d"  c(Wx  ^y  '^yy) 


i+(Sl)»  *» 

UZ 

1 

V  sin  (2  tan-1  (g))  ) 

l  2 cos  (tan-1(g))  (-sin  (tan-1^g 


V 


“ts  w!1 


\ 


2cos(tan  Kg))  (- sin  (tan  KSv,1+(^} 
cos2  (tan-1(g 


i+(f  )2  «)2 

1  1  ~B  (d.3) 


JL_  ijiH 

^E.12  «s  w 


) 


Defining  the  corresponding  nominal  angle  8n(t,  x )  =  tan-1(g^’%  j),  and  re-arranging  the 
terms  of  Eq.  (B.3),  we  get  the  linearized  version  of  the  geometric  heat  -equation  around  a 
nominal  value: 

^  «  AGllElinu(t,  x  )  =  sin 2(8n(t,  x ))  uxx(t,  x )  -  sin(20n(f,  x ))  uxy(t,  x  )  + 
at 

COS  2(9n(t,x))  Uyy(t,  X  )  +  c(-Uy(t,x)  ux(t,x)  +  u*(t,x) 


u. 


where  c 


(«S)2+(«g)*  [sin(2r)«T  -  uyy)  -  cos(20n)2<J . 
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1 


Figure  B.l:  Points  of  the  zero-level  set,  i.e.  initial  contour  [X (t),y(t)),  at  time  t,  is  shown 
on  the  left.  Those  points  whose  sample  realizations  result  in  an  average  value  of  zero  at 
time  t  +  5t  (u(t  +  5t,x)  =  Ex  [u(X  (t))]  =  0)  form  the  new  contour(T ( t  +  5t,),y(t  +  5t)) 
(on  the  right). 


T 


Figure  B.2:  Brownian  Motion  on  the  tangent  direction,  and  corresponding  interpolation 
on  square  grid. 


Figure  B.3:  Equivalent  random  walk  on  the  tangent  direction  implemented  on  the  level 
set  function  u(x,y).  The  tangent  direction  is  estimated  directly  from  the  level  set  set 

function  :  6t  =  tan-1  The  level  set  function  is  on  a  250  x  250  grid,  5t  =  0.25. 
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Figure  B.4:  Middle  row:  Brownian  Motion  on  T  is  showifltto  produce  similar  results 
with  those  in  Bottom  Row:  Geometric  Heat  flow.  The  speeds  of  the  two  algorithms  are 
different.  The  level  set  function  is  on  a  191  x  221  grid,  5t  =  0.25. 


□  □ 


□ 


Figure  B. 5:  (a)  Initial  set  of  shapes  (b)  Flow  Ct  =  cos2(2 9)kN  (c)  Flow  Ct  —  sin2(20)/cIV  . 


Figure  B.6:  (a)  Initial  set  of  shapes  (b)Flow  Ct  =  cos2(3 6)kN  ,  which  tends  to  produce 
hexagons,  (c)Flow  Ct  —  sin2(1.5(0— 7r/2 ))kN  ,  which  tends  to  produce  triangle-like  shapes. 
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Actual  Shape 


t=40 


t=80 


t=160 


Figure  B.8:  Each  row  corresponds  to  a  curve  evolution  method  with  different  n,  1st  row 
Ct  =  kN  ,  2nd  row:  Ct  =  cos2(29)kN  ,  3rd  row:  Ct  =  sin2(49)nN  . 
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Figure  B.10:  (a)  Clean  building  image  (b)  noisy  building  image  (c)  Geometric  heat  flow 
ut  =  u#  (left-right)  t  =  10,20,40  (d)  Flow  ut  =  cos2(2 9)  u#  (left-right)  t  =  10,20,40. 


Figure  B.ll:  (Top)  Diamonds  image  (Middle  row)  Geometric  heat  flow  ut  =  (Bottom 
row)  Flow  Ut  =  sin2  (26)  u#. 
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Figure  B.12:  (Top)  An  image  from  mars  pathfinder,  (First  column)  Geometric  heat  flow 
ut  =  iitf,  (Second  column)  Flow  ut  =  cos2(30)  u 


image;  2nd  row:  geometric  heat  flow 
u  =  cos2(20)ufe,  (left  to  right)  t  =  20/ 
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Figure  B.14:  Top:  House;  Left:  Geometric  heat  flow  ut  =  u^,  (top  to  bottom)  t 
40,80, 160;  Right:  Flow  ut  =  cos2(20)  u#,  (top  to  bottom)  t  =  40,80, 160. 
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